Robust optimization for casualty scheduling considering injury deterioration and point-edge mixed failures in early stage of post-earthquake relief

Objective Scientifically organizing emergency rescue activities to reduce mortality in the early stage of earthquakes. Methods A robust casualty scheduling problem to reduce the total expected death probability of the casualties is studied by considering scenarios of disrupted medical points and routes. The problem is described as a 0-1 mixed integer nonlinear programming model. An improved particle swarm optimization (PSO) algorithm is introduced to solve the model. A case study of the Lushan earthquake in China is conducted to verify the feasibility and effectiveness of the model and algorithm. Results The results show that the proposed PSO algorithm is superior to the compared genetic algorithm, immune optimization algorithm, and differential evolution algorithm. The optimization results are still robust and reliable even if some medical points fail and routes are disrupted in affected areas when considering point-edge mixed failure scenarios. Conclusion Decision makers can balance casualty treatment and system reliability based on the degree of risk preference considering the uncertainty of casualties, to achieve the optimal casualty scheduling effect.


. Introduction
Earthquakes often occur randomly without warning and bring devasting damages to affected areas. Earthquakes not only cause serious economic losses but also lead to many people's injuries or even death (1,2). Statistically, the 1976 Tangshan earthquake caused 242,769 deaths and 435,556 injuries. The 2008 Wenchuan earthquake killed 69,227 people, injured 374,643 people, and left 17,923 people missing. After a large-scale destructive earthquake, the number of casualties increases rapidly, the injury states are complicated, and the affected areas are wide. The scientific and effective treatment of the casualties is the primary task of the post-disaster relief work, and its efficiency is related to the success or failure of the post-disaster rescue. Especially the early stage of post-earthquake relief is the critical period for emergency rescue and disaster relief (3). How to organize rescue activities plays a crucial role during this golden rescue period, such as transporting and treating the casualties (4). Consequently, utilizing existing medical resources to build an efficient and reliable emergency logistics network to transport the casualties to medical points efficiently for treatment is an urgent problem to be resolved in the early post-earthquake period.
The problem of casualty transportation scheduling was usually described as an emergency facility location-allocation problem or a delivery vehicle routing problem (5)(6)(7)(8). On this basis, some studies considered the type of casualties (9,10) or the deterioration of the wounded (11) to further optimize the casualty scheduling. However, the above literature generally contained an implicit assumption that the facility was completely reliable without disruption. This assumption did not correspond to reality. Especially in the post-earthquake emergency logistics system, the facilities and routings are easy to fail due to the impact of earthquake disasters and their secondary disasters. The disruption factors were considered in some literature (12)(13)(14).
To sum up, researchers have carried out in-depth research on the optimization of emergency scheduling for casualty transportation after earthquake disasters, but there are still some problems that are required further study. Firstly, most of the existing research studied the problem of injury deterioration and casualty scheduling alone without considering the two problems simultaneously. Secondly, most of the current studies assumed that the road and facilities in the earthquake disaster logistics network were completely reliable without failures. Although some literature considered road disruption or facility failures, there are few works considering the mixed failure state of points (facility) and edges (route) in the emergency logistics network simultaneously. Thirdly, existing studies usually assumed that information such as the number of casualties was determined or described as fuzzy parameters. In reality, such a number cannot be estimated exactly. To fill up this gap, this paper considers the point-edge mixed failure scenario, combines the evolution of injury situation, adopts the robust optimization method to deal with the uncertain amounts of casualties, and further studies the optimization problem of casualty transport scheduling in the early stage of post-earthquake relief.
The main differences between this paper and previous studies are summarized in Table 1. More specifically, the contributions of this paper are as follows. (1) We propose a new casualty scheduling robust optimization problem for post-earthquake relief. As mentioned in the previous subsection, the proposed problem considers some new characteristics in the early stage of postearthquake relief, such as point-edge mixed failure scenarios, injury deterioration, and multiple transportation modes. (2) A mixed integer non-linear programming (MINP) model is established to formulate the proposed problem. (3) According to the characteristics of the model, an improved integer-encoded particle swarm optimization (PSO) algorithm is designed in this paper. The improved PSO is compared with the genetic algorithm (GA), immune algorithm (IA), and differential evolution algorithm (DE).

. Literature review
This section will review relevant literature in the following three categories.

. . Traditional casualty scheduling problems
The problem of casualty transportation scheduling is usually described as an emergency facility location-allocation problem or a delivery vehicle routing problem (5)(6)(7)(8). Mansoori et al. (15) proposed a multi-objective humanitarian supply chain design problem that minimizes the total number of injured not transferred to hospitals. Fiedrich et al. (26) proposed a dynamic optimization model where the total number of fatalities during the initial search-and-rescue period after strong earthquakes is minimized. Andersson et al. (27) established a support tool for dynamic ambulance relocation and automatic ambulance dispatching. Xie et al. (28) formulated a lane-based evacuation network optimization problem that integrates lane reversal and crossing elimination strategies. A Lagrangian relaxation algorithm based on the principles of Tabu search is designed to solve the model. Mclay et al. (29) proposed an improved Markov decision model to optimize ambulance dispatch dynamically to maximize the number of critically injured patients. Knyazkov et al. (30) studied the present situation of emergency transport for patients with acute coronary heart disease in St. Petersburg. The optimization of ambulance routes was proposed according to the road network in the city, to improve the number of ambulances per capita and the speed of ambulance response in this work. Sung et al. (31) transformed the problem of treating the injured after a disaster into the dispatching problem of emergency ambulances with the goal of maximizing the expected survival rate and designed a column generation algorithm. Repoussis et al. (32) proposed a mixed integer programming formulation for the combined ambulance dispatching, patient-to-hospital assignment, and treatment ordering problem. The objectives are to minimize the overall response time and the total flow time required to treat all patients. Shavarani et al. (33) described how to properly allocate existing emergency vehicles to hospitals and effectively plan vehicle routes after a disaster in a densely populated area, to maximize patient survival. Sun et al. (16) proposed an emergency model of the location-transportation-allocation problem. The objective is to minimize the total cost and the sum of injury severity scores. The integrated research on the problem of casualty scheduling and emergency facility location is studied in some literature (34,35). Sheu et al. (36) proposed a method for designing a seamless centralized emergency supply network. A three-stage multi-objective (travel distance minimization, operational cost minimization, and psychological cost minimization) mixed-integer linear programming model is built to describe the problem. Hu et al. (37) proposed a mixed integer programming model that considers the uncertainty of the number of injured people and integrates the decision of locating shelters and transferring injured people. Chou et al. (17) proposed a patient transportation and assignment model considering the routing of ambulances and operational conditions of hospitals.
. . Casualty scheduling problems considering wounded types A large number of wounded personnel emerged in the early stage after an earthquake, thus it is important to dispatch the wounded personnel considering the wounded type (9, 10).

. . Casualty scheduling problems considering disruptions
The facilities (points) and the roads (edges) of an emergency logistics network might be disrupted after a large-scale earthquake. The logistics network design problem considering facility failures is studied in some literature. Bayram et al. (45) studied the shelter location and evacuation route assignment problem considering the disruption/degradation of the evacuation road network structure. Cheng et al. (23) adopted a two-stage robust optimization framework to study the robust fixed cost location problem in the case of uncertain demand and facility disruption, and developed a column constraint generation algorithm to accurately solve the model. Zhou et al. (46) constructed a location-allocation model for emergency facilities suitable for the initial stage of post-earthquake rescue, considering facility disruption and multiple types of fuzzy demand. Mohammadi et al. (47) studied a multi-objective reliable optimization model to organize a humanitarian relief chain. They made a broad range of decisions, including reliable facility locationallocation, fair distribution of relief items, assignment of victims, and routing of trucks. Sun et al. (24) proposed a scenario-based robust dual-objective optimization model to study the location of medical facilities, casualty transport, and relief material distribution under the temporary medical point failure scenarios. There are also some studies on the design of emergency logistics networks considering edge (road network) failures. Sabouhi et al. (48) proposed a comprehensive stochastic programming model for the distribution of relief materials in disaster areas, taking the demand

. Description of injury deterioration
The injury deterioration can be described by a Markov chain proposed by Wilson et al. (43). The injury state of the wounded evolves toward the direction of gradual aggravation or continues to maintain the original state, with a certain probability in the process of waiting for rescue. And death is the end of the injury evolution in process. If the injury at the current stage is severe, the injury may be transferred to the death state or remain severe if effective treatment is not available at the next stage. If the injury is minor at the current stage, the injury may be transferred to severe or remain minor if no rescue is available at the next stage. It can be seen from the evolution of the injury state that the state transfer is stochastic, and the state of the wounded at the next stage only depends on the current state independent of the historical state, which is consistent with the no aftereffect of the Markov process. Therefore, the evolution of the injury has a Markov character.
In this paper, the evolution process of injury is divided into three finite-state Markov processes: minor injury (Green, G), severe injury (Red, R), and death (D). The process is repeated until the casualty dies or is saved. Let the initial time t = 0. All rescue teams are ready at the initial time. The injury level of the casualty is randomly and unidirectionally changed once per minute without any medical treatment, as shown in Figure 1.
It is assumed that there is a linear relationship between the severity of injury and the time before the wounded is treated. According to the random transfer probability between the injury states, the possible death probability of the casualty j whose initial injury state is R is expressed as formula (1).
In expression (1), T r ij indicates the time required to transport a severe casualty from the affected point i to the medical point j. L r represents the injury state of the casualty in the initial time is r. p rd indicates the probability that a casualty's injury state changes from severe to death and let p rd = 0.1 in this article referring to literature (43).
If p rd T r ij > 1, the wounded may die in a waiting process. The death probability of the severe casualty while waiting for rescue can be expressed by Equation (2).
The death probability of the minor casualty can be written by the transition probability equation. That is, the probability of changing from a minor injury to a severe injury is calculated first, and then the probability of deteriorating from the severe injury to death is calculated. The time T that the wounded fully evolves from G to R should be first judged when calculating the death probability. According to Figure 1, the transition probability from state G to R is 0.05. Let 0.05 * T = 1, then, T = 20 when G evolves into R. When the ambulance arrival time is <20 units, the wounded does not evolve to death, and the death probability is 0. When the arrival time of the ambulance is >20 units, the death probability in this . So, the probability function of death in a minor state is a time-segment function, see expression (3).
Therefore, the death probability functional of a casualty can be expressed as (4).

. Problem description
The impact of aftershocks, mudslides, and other secondary disasters, not only leads to medical points failure but also causes road damage or interruption, in the early stage of post-earthquake relief. There is a mixed failure of facility nodes and routes in the emergency logistics network. The transport network diagram with the point-edge mixed failures in the early stage of post-earthquake relief is shown in Figure 2. As can be seen from Figure 2, based on the failures of facility points H and G, the secondary disaster also leads to the interruption of roads between affected point 4 and medical point D, and between affected point 4 and medical point E. At this time, ambulances cannot normally pass on the two roads, so it is necessary to use helicopters to transport the wounded to ensure timely treatment, reduce the death rate of the wounded, and improve the efficiency of emergency rescue.
The purpose of emergency rescue is to treat as many casualties as possible in the shortest time and reduce the death rate of the wounded. Therefore, the goal of the model is to minimize the expected death probability of casualties.

. . Assumptions
Model assumptions are as follows. First, the location of affected points and medical points is known, and the resources available at the medical point are limited. Second, each medical point can serve multiple casualty groups at the same time, and each casualty group can be assigned just one medical point. Third, the casualty groups have been classified in advance according to the trauma index of the injured. The severe injuries should be treated first according to the medical resources and ambulance and helicopter transport conditions. Four, each ambulance can only transport two casualties at a time, and each helicopter can transport four casualties at a time.

. . Notations
The parameters are as follows. ε ij : Degree of route damage between affected point i and medical point j; if ε ij > 0.5, the route is disrupted, and the ambulances cannot pass. Then, the helicopter transportation should be considered. With the support of modern remote sensing technology and communication technology, the damage degree of route in the early stage after earthquake can be sensed in real time without constant trial and error. d ′ ij : Generalized transportation distance from affected point i to medical point j. The more damaged the route, the slower the vehicle speed, and the longer the transport time. That is, the longer the d ′ ij . d ′ ij can be obtained according to the method mentioned in literature  The decision variables are as follows. k ijw : The number of casualties with injury state w transported from affected point i to medical point j; X ij : A binary variable. If medical point j allocated to affected point i, X ij = 1;otherwise X ij = 0; X i′j : A binary variable. If medical point j reallocated to affected point i ′ after the initial allocated medical point j ′ failed, X i ′ j = 1;otherwise X i ′ j = 0; Y i ′ jh : A binary variable. If ambulance h travels from affected point i ′ to medical point j after the initial allocated medical point Y i ′ jh : A binary variable. If ambulance h travels from affected point i ′ to medical point j after the initial allocated medical point j ′ failed, Y ijh = 1; otherwise Y ijh = 0; Z ijn : A binary variable. If helicopter n travels from affected point i to medical point j, Z ijn = 1; otherwise Z ijn = 0; Z i ′ jn : A binary variable. If helicopter n travels from affected point i ′ to medical point j after the initial allocated medical point j ′ failed, Z i ′ jn =1; otherwise Z i ′ jn = 0.

. . Mathematical formulation
The mathematical model can be formulated as follows.
The objective function (5) is to minimize the total expected death probability of casualties. Constraint (6) denotes the death probability of casualties with different injury states. Constraint (7) is the resource limitation of medical points. Constraint (8) denotes the initial resources of medical points. Constraint (9) denotes the restriction of available resources at medical points at time t. Constraints (10) and (11) denote the restrictions on the number of available ambulances and helicopters from affected points to medical points respectively. Constraints (12) and (13) denote the waiting time for casualties with injury state g and r, respectively. Constraint (14) denotes the transportation distances from affected points to medical points. Constraint (15) denotes the time taken to transport casualties from affected points to medical points. Constraint (16) denotes each affected point will be served by the nearest open backup medical point if its initial allocated medical point has failed. Constraint (17) is a 0-1 variables restriction.

. . A robust optimization model
As the destructive earthquake occurred suddenly, and aftershocks may result in injuries in the initial stage after an earthquake, the number of casualties in affected points cannot be accurately estimated. Therefore, robust optimization was used to deal with the number of casualties to reduce the risk caused by uncertainty.
The parameter Ŵ iw and corresponding variables are introduced to process the objective function and constraint Referring (51). Objective function equation (5) has uncertain variables of f iw .
Define the value range off iw as f iw −f iw , f iw +f iw . f iw is the nominal value of uncertain number of casualties.f iw is the maximum deviation of the uncertain number of casualties.
A protection function for the number of casualties is introduced here as equation (20). Ŵ iw is the control coefficient of the number of casualties. A protection function for the number of casualties is introduced here as equation (20). Ŵ iw is the control coefficient of the number of casualties. Ŵ iw ∈ [0, |J iw |]. Where J iw is the number of Ŵ iw . The objective of this robust method is that the number of casualties varies within their interval at most ⌊Ŵ iw ⌋ affected points. Ŵ iw is the balance between robustness and optimality. The bigger the Ŵ iw , the more conservative the mode. The protection function for the number of casualties is introduced as with ϕ(X, Ŵ iw ) consideration of the uncertainty of the number of casualties. The objective function (5) can be transformed into formula (20).
Where, S iw represents the group set of the maximum number of casualties deviating from the nominal value. When Ŵ iw = 0, ∀i, w, the robust model is equivalent to the nominal model.
If Ŵ iw = 0, ∀i, w,all the number of casualties are nominal. When Ŵ iw = γ , the number of casualties in all injury states deviates from the nominal values, the model is equivalent to the Soyster model. As Ŵ iw changes, the conservatism of the model also changes accordingly. The objective function (5) can be finally transformed into the following expressions based on the strong duality. . Algorithm The built model is a 0-1 mixed integer non-linear programming model, which cannot be solved by exact algorithms such as branch and bound algorithm or operations research software such as CPLEX. An improved integer-coded PSO algorithm is designed to solve the model.
The algorithm steps are described as follows.

. . Population initialization
x Initialization of parameters: Set population size as popsize, maximum iterations as max gen, learning factors as c 1 and c 2 , inertia weight as w.
y Particle encoding and decoding: Supposing there are n casualty groups and k medical points in the model, each particle has a code length of n. Each position of the particle is a positive integer randomly generated between 1 and k, and denotes the relation of assign between a casualty group and a medical point. Taking Figure 3 for example, there are 6 medical points in the affected area, which are required to provide rescue services to 11 casualty groups. The length of the particle is 11. Casualty group 2 and 11 are assigned to medical point 1, casualty group 1 and 10 are assigned to medical point 2, and so on.
z Initialization of the best value of individuals and groups: The initial location and velocity of particles can be randomly  generated. The fitness value of the current population can be calculated by the fitness function. The fitness function is the objective function (25). The individual position is the optimal position of the current individual P best . The minimum value of the current particle P best_value is the initial optimal value of individuals that can be got by comparing the optimal values of all individuals. Then set the initial P best_value as the best value of groups g best .

. . Updating speeds and locations
Population speed and location can be updated as expressions (26) and (27).
In which, r 1 ,r 2 are uniform random numbers within the range of [0,1].
Traditional PSO is mainly used to solve the continuous optimization problem, while the model in this paper belongs to a discrete combinatorial problem. So, the updated particle positions and speeds need integer processing. Figure 4A is the result of an updating particle according to the particle update equations, and Figure 4B is the result of a particle rounds up to integers. Here, casualty group 5 and 10 are assigned to medical point 1, casualty group 1 and 8 are assigned to medical point 2, and so on.

. . Updating the best value of individuals
Comparing the current particle fitness value fitness with the individual historical optimal value P best_value , the smallest value is taken as the current individual optimal value P best_value , and its corresponding position is taken as the current individual optimal position P best .

. . Updating the best value of groups
Comparing the current individual optimal value P best_value with the group historical optimal value g best_value , the smallest value is taken as the current group optimal value g best_value , and its corresponding position is taken as the current group global optimal position g best .

. . Termination of the algorithm
Determining whether the algorithm reaches the maximum iteration, if so, the algorithm terminates and outputs the optimal solution. Otherwise, go to step 4.2 and continue the iteration.

. Computational experiments . . Case description and data
The data of emergency medical rescue for Ya'an, in the Lushan earthquake, Sichuan province of China is used for numerical simulation. Taking town or township as units, Lushan County has jurisdiction over 5 towns and 4 townships. So, 9 affected points are set in this article. The nominal number of casualties in each affected point is shown in Supplementary Table S1. The number of treatment resources required by casualties with injury state r and g is 3 units and 2 units respectively.
Data shows that the Chinese government and military deployed 15 helicopters to transport and rescue casualties in the Lushan earthquake. Therefore, the number of available helicopters in medical points is set as 15. According to the basic Standard for Medical Institutions (Trial) issued by the Chinese Ministry of Health, there is one ambulance for every 50,000 people in the city. The number of ambulances in Ya'an city and Chengdu city is about 30 and 300, respectively. In this paper, it is assumed that 50% of ambulances in Ya'an city can be used normally, and due to the time crunch, only 50% of ambulances from nearby Chengdu city are used for the transport of casualties. So, the total number of available ambulances in this paper is set as 165. The other parameters to the model are set as following. C r = 3 units, C g = 2 units, v h = 40 km/h, v n = 120 km/h, q h = 2 persons, q n = 4 persons, T = 30 h. The number of medical resource limitations for each medical point is shown in Supplementary Table S2. The distance between medical points and affected points is shown in Supplementary Table S3.   Referring to literature (46), four failure scenarios are set in Supplementary Table S4

. . Calculation results
The PSO parameters are set as follows. Population size popsize = 100, learning factors c 1 = c 2 = 2, inertia weight w= 1, the maximum iterations max gen = 300. MATLAB R2018b is used for programming, and it runs on a notebook computer with Intel(R) Core (TM) I5-5200U CPU and 12G memory. The optimal solution of the results of 10 operations is taken as the final solution. The convergence time of the algorithm is 10.171 s, and the objective function is 531.629. The convergence diagram of the algorithm is shown in Figure 5. Table 2 and Supplementary Table S6 are treatment information at medical points and affected points respectively.

. . Considering failures vs. without considering failures
To test the reliability of the scheme got by considering the point-edge mixed failures, a comparison was made between two schemes. Scheme 1 is the current optimization results got by   (Table 3).
The results show that when there is no failure (scenario 1), the casualty scheduling results under scheme 1 are slightly worse than those under scheme 2. However, when the failure occurs after the earthquake, the solution result considering the failures is better than that without considering the failures. As shown in scenario 3, the total expected death probability of casualties under scheme 2 was 60.22% higher than scheme 1. Therefore, considering the point-edge mixed failure in advance can make more casualties receive timely treatment and reduce the death rate. The reliability of the system can be improved by considering the failures in the design stage of the emergency rescue network.

. . Robust optimization vs. deterministic optimization
To verify the effectiveness of the robust optimization model, the running results of the robust optimization model and the deterministic model were compared in the same scenario. Four numerical examples are designed according to different control coefficients Ŵ w , and the maximum disturbance value of casualty number deviating from the nominal value is 20%, that isf iw = 0.2 × f iw . Table 4 shows the sensitivity analysis results of different robust control coefficients Ŵ w . In general, although the objective function value of the robust optimization model is higher than that of the deterministic model, the gap is not significant. It shows that the robust optimization model can reduce uncertain risk. From the results of sensitivity analysis of Ŵ w , the smaller the Ŵ w , the stronger the robustness of the model. With the increase of the uncertainty of casualty numbers, the objective function value increases. The gap between the robust optimization model and the deterministic model (when Ŵ w is 0) is gradually increasing. The appropriate scheme should be chosen according to the uncertain situation in practical decision-making.

. . Algorithm performance analysis
To test the performance of PSO, the algorithm was compared with GA, IA and DE. Each algorithm was run 10 times. The results are shown in Table 5 that the PSO is significantly superior to other algorithms. The average objective function value calculated by PSO reduces by 53.2% compared with GA, 37.5% compared with IA, and 53.0% compared with DE. The computational accuracy of PSO is better than GA, IA and DE, as shown in Figure 6.

. Conclusions
A robust optimization model was established in this paper to minimize the total expected death probability of casualties in the Frontiers in Public Health frontiersin.org . /fpubh. . early stage of emergency rescue after an earthquake, considering the characteristics such as casualty classification, injury deterioration, point-edge mixed failures, medical resource limitation, and casualty number uncertainty. An improved PSO is proposed to solve the problem. Numerical experiments are conducted to verify the effectiveness of the model and algorithm under the context of the Lushan earthquake in China.
The results reveal that the operation effect of the emergency rescue network is significantly improved, and the optimization results are more reliable if the failure scenario is considered in advance. Moreover, robust optimization considering the casualty uncertainty can reduce the uncertainty risk of the system. Therefore, it is necessary to consider the point-edge mixed failures and casualty uncertainty simultaneously in the design stage of the emergency rescue network to establish a more reliable emergency rescue network.
Future research can consider the characteristics of the later stage of post-earthquake rescue comprehensively to study the casualty scheduling problem in the post-earthquake emergency recovery period while considering factors such as resuming normal passage of affected points and transporting the injured.

Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.